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ylBSTRACT 


Recursive least-squares estimates for processes that 
can be generated from finite-dimensional linear systems are 
usually obtained via an nxn matrix Riccati differential 
equation, where n is the dimension of the state space. In 
this new recursive method the gain matrix for the Kalman 
filter and the convariance of the state vector are computed 
not via the Riccati equation, but from certain other equat- 
ions, These differential equations are said to be of 
Chandrasekhar- type , because they are similar to certain 
equations introduced in 1948 by the astrophysicist, 

S. Chandrasekhar, to solve finite-interval Wiener-Hopt 
equations arising in radiative transfer, Chandrasekhar 
extended Ambarzamian’ s invariance principles to solve the 
above problem and in fact in recent years this served as a 
stimulas for much activity by Bellman in the development of an 
invariant approach to the solution of various transport process. 
The "invariant imbedding" idea resulted in the reduction of 
the basic boundary value problem of transport theory to an 
equivalent initial value system, a significant computational 
advance , 

Initial experience has shown that there is some 
computational savings in the new method and the loss of 
positive definiteness of the covariance matrix is less vulner- 
able. 
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CHAPTER I 


INTRODUCTION 

The work of R. E. Kalman is a most significant contribution to 
filtering and prediction theory, since the original work of Weiner 
(8). The Kalman filtering theory introduces a new look at the 
classical problems of prediction, smootfiing, and filtering. More 
specifically, Kalman's method has the following features: 

1. The linear dynamic system is described by the state variables and 
state equations. This not only represeists a modern approach to the 
systems problem, but also makes machine computation simpler. 

2. The Kalman filtering theory treats stationary and nonstationary 
random processes, single-variable and multivariable systems, all in a 
unified manner. 

The Kalman filtering problem can be stated in general as: given 
y(t) = Z(t) + v{t), where y(t) is a mebsage in the form of a signal 
corrupted by additive noise and Z(t) is the actual signal and v{t) 
is the noise, determine the value of the Z(t) in the sense of the 
minimum mean-square error. The error is defined to be the difference 
between the actual output of the filter and the signal component of 
the input message. 

Generally to determine the value of Z(t) at sometime t = tj, 
given the measured or observed value {y{T), t^ < t < the time 

tj can be less than, equal to, or greater than tj^. These three cases 


can be defined as 
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tj < \ 

smoothing (interpolation) 

11 

filtering 

> tk 

predicting 


The block diagram illustrating the general philosophy of the 
Kalman filter is shown in Fig. 1. 



Kalman 



(x > 0) Smoothing (Interpolation) 
M^j(s) = 1 Filtering 

M^i(s) = (t > 0) Predicting 


FIG. 1 KALMAN FILTER 
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This research will include the study and comparison of a new 
algorithm for recursive state estimation via Chandrasekhar-type 
equations (R). 

By this method, the gain matrix for the Kalman filter is 
obtained directly, without having to solve separately for the error- 
covariance matrix. In general* the gain matrix is obtained by solu- 
tion of n(n+l)/2 simultaneous nonlinear differential equation of 
Riccati-type, where n is the dimension of the state space. But in 
this new method, it only requires the solution of n{m+p) simultaneous 
nonlinear differential equations, where m and p are the dimension of 
the input and observation processes, respectively. In most practical 
cases n> p + m, and our experience shows that whenever n> p + m, there 
Is some computan onal saving in using this new method. 
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Notation 

The “dot" notation will be used for derivatives (i.e., ^ = x); 
and "prime” will be used for transpose of a matrix as shown (F'). 
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CHAPTER II 

F0RI4ULATI0N OF PROBLEM AND SOLUTION TECHNIQUE 
Problem Statement 

Consider the standard Kalman state-space model where the problem 
Is to calculate linear least-square estimates of a signal process 
2 (*J from p- vector observations of the form 

y(t) = Z't) + v(t) ttto (2-1) 

where 

Elv(t)v'(s)} = Ip fi(t - s) (2-2) 

and z{’) Is given by a state-space model 

Z(t) ^ HX{t). t > to (2-3) 

X(t) = FX(t) + GU(t). X(to) = Xq (2-4) 

»:^here x(-) is an n-vector, U(*)is a m-vector, and 

EIXq] = 0, EIXqXJ] = Hq. EtU(t) X^l = 0 (2-5) 

the input process u(-) is white with covariance function 



ElU(t)U'{sr = Q6(t - s). ElU(t) v'(s)] = C6(t - s) (2- 6) 

The matrices F, G, H» Q, C are assumed to be known and for simplic- 
ity we also assume that 

H has fiill rank p (2-7) 




] 


It is desired to calcula+j 
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A 

x(t) “ the linear least-square estimate of x(t), given {y (t), 

t(j < T < t). 

t ' 

By linearity it follows 'Jiat 

z(t) = HX(t) (2-8) 

The Kalman filter for this problem 1s 

X(t) = FX(t) + K(t)(y(t) - HX(t)), X(t„) = 0 (2-9) 

where 

K(t) = P(t)H" + GC (2-10) 

The form of this best linear estimator, which is called the 
Kalman filter is shown in Fig. 2. (7) 

The nxn matrix P(*) is the mean-square error in the estimated of 
the states 

P(t) = E[x(t)-x(t)][x(t)-^(t)]" (2-11) 

This matrix is computed as the unique solution of the non-linear 
matrix differential equation of Riccati type. 

Rt) = FP(t) + P(t)F' - K(t)K"{t) + GQg". 

P(tj,) = Ho (2-12) 

since P(*) is syumetric, this Riccati equation involves the solution 
of n(n+l)/2 simultaneous nonlinear differential equations, which has 
generally to be carried out on an analog or digital computer. The 
fact, clear fron (2-11), that P(*) has to be nonnegative-definite 
often places certain strigent accuracy requirements on the computation 
and may require special attention. However, by this method K(*) is 


/ 










calculated not via Riccati equation but from certain equations, which 
are called of Chandrasekhar- type. 

Chandrasekhar-Type Equations 

The Chandrasekhar-type equations are as follow (5): 

Ut) = .[Yi(t)Y^(t) - Y2(t)Y2(t)]H\ K(to) = HoH"+ GC (2-13) 

Y^(t) = [F - K{t)H}Y-|(t), Yi(to) - i = 1»2 (2-14) 

The initial condition matrices are to be determined from the 
following procedure 

d^fhq + ^ 0 ^' + gqg"- (HoH' + gc)(tqH" + gc)' ( 2 - 15 ) 

Let us assume that D has rank o, and a < n. From (2-15) i:: is 
clear that D is a symmetric matrix and it will have real eigenvalues. 
Here in various ways, which will be described later, we can write D 
as 

D = [LiLJ - I 2 L 2 I (2-16) 

If we assume that 

$ = the number of positive eigenvalues of D. 

Then the dimensions of matrices L, would be 

Li = nxe, L 2 = nx(a - p) (2-17) 

where the number of positive eigenvalues of D is to be calculated, 
however, there is no need to calculate the eigenvalues. 

The best way to find the initial conditions is to use 
Cholesky Decomposition which will be described later in this chapter. 


/ 


10 


Therefore let us make a so-called lower tn’angular-di agonal -upper 
triangle (LDU) decomposition of D as 

D« LcSLq = [Lj L2 ]S[Li L^] (2-18) 

where we define S to be the a x a signature matrix of 0 
S “ diagd ,1 . 1 . . . ,1 ,-l . ,-l } 

Now by these assumptions and with (2-18) we can write the 
Chandrasekhar- type equations more compactly as 

iC(t) = Y(t)SY"(t)H", K(to) = KqH" + GC (2-19) 

Y(t) = [F - K(t)HlY(t), YCtjj) = [Lj L^] (2-20) 

where 

K(*) is an (n x p) matrix 

Y(*} is an (n x ct) matrix 

F is an (n x n) matrix 

H is an (p x n) matrix 

The important aspect of these two last equations (2-19), (2-20) 
is that the Kalman gain K(-) can be found directly without going 
through the matrix ?(*)• Note that if it is desired to find P(*)» we 
may use (2-10) to write 

p(t) = Y(t)SY'(t), P(to) = Ho (2-21) 

We would like to develop the new algorithm based on the 
Chandrasekhar- type equations (2-13) and (2-14). Here we have n(a + p) 
simultaneous nonlinear equations. In many practical cases a is less 
than or equal to p or to m, where mis the number of inputs; and often 
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m » p = 1. 

Comparing n(n+l)/2 simultaneous nonlinear differential equations 
of Riccati-type with n(a+p) simultaneous nonlinear differential equa- 
tions of Chandrasekhar- type, we see that whenever n ^ p + a, we will 
have a substantial computational saving in using the Chandrasekhar- 
type equations (2-13) and 2-14). Good experience wh^ch has been 
obtained in several examples, shows that the equations (2-13) and 
(2-14) are numerically well behaved. 

Let us now define, the steady state behavior (t •> «) of K(*) as 

follow 

K=limK{t) (2-22) 

t ">■ » 

Because (2-13) and (2-14) are a set of simultaneous equations, 
it is computationally preferable to find K by computing the solution 
of Chandrasekhar- type equation until the solution remains fairly 
constant. This will be shown on an *'*xampTa later in this chapter. 
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Determination of the Initial Condition Matrices Li and L 2 

CKOLESKY DECOMPOSITION 

Any symmetric positive semidefinite n x n matrix D may be writ- 
ten in the factored form (6): 


“ ■ 


* « 

• • 

Dll 


Lii 0 ... 0 

Iji L21 ... L^i 

Di 2 D22 

« « 


C21 L22 0 

• • • 

0 L22 

■ f • 4 

• « 

« m 


• 4 » 

• ■ 4 

0 

4 4 4 4 

• 4 4 

t^in ... ^nn 

^ m 


*-ni . . ,W 

• ^ 

0 0 ... Lpp 


This is called lower triangle-upper triangular decomposition, 
Cholesky gave the following recursive algorithm for computing L, 
.For i = 1, n. 


i-i 2 




13 


Example: 1 



Exampl e : 2 


[3 2 12j 


where o = 3, e = 2 
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*1 2 si Ln 0 0 1 ^ ® 

2 8 2 = L21 L22 0 0 1 0 0 L22 L32 

3 2 12 L31 L32 L33 _0 0 ° 

. ** 

lifter multiplying the right-hand side and equating the corres 
ponding elements 



Form Column 1 of L 
Lii =jDn/Si 
L 21 = D 21 /S 1 L 11 


L31 ^ D31/S1L11 



lZ2 - 
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The fortran algorithm is as follows: 
L(l,1) - SQRT(ABS(D(l.l)/s(l))) 

DO 100 I=2,N 

L(MH 

IF{L(1,1).EQ.0. )L(I,1)=0 
,TF(L(M,EQ,0.) GO TO 100 
L(I,l)=D{I.l)/(sn)*Ul,l)) 

100 CONTINUE 

DO 500 0=2. N 

JM1=J-1 
0P1=0+1 
L(J,0)=D(J.J) 

DO 200 K=1 , OMl 

200 L(J,J)=(L(J,J)-S(K)*L(J.K) ** 2 ) 

L(J.J)=SQRT(ABS{L(J.0)/S(0))) 
IF(J.EQ.N) GO TO 500 
DO 400 r=OPl,N 
L( 0 ,I)= 0 . 

IF(L(0,J).EQ.0.)L(J,J)=0 
IF(L(J.0).EQ.O)L{J,O)=O. 
IFa(J.J).EQ.0)60 TO 400 
Ui.oHCi.o) 

DO 300 K=1 , JMl 
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300 LCI,J)=(L(I..^)**S(K)*L(J.K)*U(I.L})/S(J) 

L(I,J)=L(I,u)/L(J,J) 

400 CONTINUE 
500 CONTINUE 

We now turn our attention to investigate the variety of special 
cases. These special cases are of great important. 

Case I : (Low Initial Uncertainty) 

Let us now assume that the initial state x(to) is fairly known. 


that is we may write 

EIXoXol = Ho 0 (2-23) 

Substitute no - 0 into the equation (2-15), then the matrix D 
simplifies to 

D = G(Q -CC")G'> 0 (2-24) 

Here D can have no negative eigenvalues, so that D may be written as 

D=LjLJ, L = nxa (2-25) 

The Chandrasekhar- type equations now can be simplified to the fewer 
than n(m+p) equations 

k(t) = Yi(t)Yj(t)H', K(t„) = GC (2-26) 

il(t) = [F - K(t)HlYj(t), Yi(to) = L) (2-27) 


If it is desired to find P(*}. the mean-square error in the estimate 
of the states, it can be found via 

P(t) = Y,(t)Y;(t); P(t,) = 0 


(2-28) 
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The above equation can also be written as 

P(t) = Y (T)Y:(T)dT ■ (2-29) 

0 ^ ^ 

A significant aspect of this last formula is that no matter how 
Inaccurately Yi(*) may have been computed, the product Yj(OVi(*) is 
always nonnegative-definite. This property may be last in the subse- 
quent quadrature (2-29), but we realize that it is much easier to be 
careful in a simple quadrature than in the solution of the Rlccati- 
type differential equation (2-12). 

Let us now consider the following example where the Kalman gain 
K(') is find via Chandrasekhar- type equations. 
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Example: 3 

Consider the Chandrasekhar-type equations 


K(t) = Yi(t) YjCt) h'. 


K(t ) ^ 6C 


Yj(t) = IF - K(t) H] Yj(t) Yj(t ) = Lj 

We would like to plot the values of K(*) - P(*)H R"^ for a 
five-state system with no == 0 and 


F ^ 


6 ^ 


-0.0297 

0.331 

-1.13 

0 

0 

f -0.0297 
0.331 
-1.13 
0 
0 


H « 


-1 

-0.0042 

0.128 

0 

1 

0 

0.381 

0.067 

0 

0 


To 0 0 1 6 

[o 0 0 0 1 


0 

•0.0461 

*0.803 

1 

0 

0 

0,040 
1.59 
0 
0 


0.0438 

0 

0 

0 

0 


Q = 


0.01 

0 

0 


0 

0.001 

0 


0 

0 

0 

0 

0 


0 
0 

O.OOl 
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and with E{v(t)v"(s)} « R 6(t - s) where 


And 


R = 


C 


*0.001 

0 

0 

0.001 

4 


’o 

0 


' 0.003 

0 

0 

\ 

-0.035 

_0 

0 

m 


0.124 


0,00 

0.00 


Substitute the values into Chandrasekhar- type equations yield: 



Kii 

■ 

Ki 2 


Y11Y41 

Y11Y51 


K21 

K22 


Y21Y41 

Y21Y51 

0.001 

ki 

K32 


Y31Y41 

Y31Y51 


K41 

i ^2 


Y41YU 

Y41Y51 



ks 2 _ 


YsiYm 

% 

Y51Y51 

1 


fi.: 


m 

-0.0297 



0.0438-Kn 

-Ki2 

Yn 



0.331 

-0.0042 

-0.0461 

“^21 

-K 22 

Y21 

Y31 

= 

-1.13 

0.128 

-0.803 , 

-K31 

“K 32 

Y31 

Y41 


0 

0 

1 

-K41 

-K 42 

Ym 

Ysi 


0 

L 

1 

0 

“K51 

-K 52 

Y51 

0 
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Where we would finally have fifteen simultaneous differential 
equations as follows: 

kn ^ 1000 Y 11 Y 41 
= 1000 Y11Y51 

K21 = 1000 Y21Y41 
K22 = 1000 Y21Y51 
K31 - 1000 Y31Y41 
K 32 = 1000 Y31Y51 
. K41 = 1000 YmYm 

K42 “ 1000 Y^jYjj 
K51 - 1000 Yg^Y^j 
K52 “ 1000 YgjYgj 

Yji “ 0.0297 Yji - Y 21 + 0.0438 Y^j - ^llYui " ^^laYsi 
Y 21 = 0.331 Yu - 0.0042 Y 21 - 0.0461 Ygj - K^jY^j - ^^ 2^51 
Y 31 = 1.13 Yu + 0,128 Y 21 - 0,803 Ygj Kg^Y^^ - K32Ygj 

Ym ” Y 3 J - K 41 Y 4 J - ^42^51 
Y5I “ Y2I - K51Y4I - Kg 2 Y 51 

I have used IBM-1130 t'.* calculate the values of Kalman gain. 

Fig. 3 shows the plot of these values. We see that the solution of 
the Chandrasekhar-type equations run quire smoothly into steady state. 

The complete computer program and numerical results, using the 
4th order Runge-Kutta method, are shown in Appendix I. 
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EXAMPLE: 4 

As a further example consider an RLC circuit as shown in Fig. 4. 
The voltage 2{t) is the out put of a system with input, the white 
noise U(t). By Kirchhoff's law the following voltage equation may be 
written for the circuit 

U(t) = L + Ri (t) + 1 idt (2-30) 

ot C 

where R, L, and C are the resistance, inductance, and capacitance, 
respectively, and i(t) is the circuit current. 

Let the variable q be the electric charge then: 

(Z-31) 

q(tj) = 0 (2-3Z) 

Equation (2-3D) may be written 

U(t) = + (2-33) 

Equation (2-33) Is a second-order equation governing the behavior of 
charge densi ty. 

To put Equation (2-33) in the standard state-variable form use is 
made of the definition of Equation (2-31). Thus Equation (2-33) is 
equivalent to the two coupled first-order equations. 


- Ri(t) - jq + u(t) 


(2-34) 
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Defining 

X, 5 q Xj s 1(t) 
Equation (2-34) may be written as 

x{t) = F X(t) + GU(t) 

where 



•r * 

h 


0 

r 

I* 

X(t) = 

X2 

. U(t) = (Uil . F = 

1 

;cc 

L 

« 

G « 


Similarly output voltage Z{t) may be written as 

Z(t) = 1 q = ^ Xj 
or 

Z(t) = H X(t) , K = [ ^ 0] , 

Suppose that our measurement is y(t) 

y(t) = 2(t) + v(t) , v(t) s observable noise 
Let us assume that the input noise and the observable noise 
uncorrelated, so that 

E [U(t) v'(S)] = 0 
and 

E[U(t) U'(s}] = [1] 


(2-35) 


(2-36) 


(2-37) 


(2-38) 

are 

(2-39a) 


Ho = 0 


(2-3%) 
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With these assumptions we can write Chandrasekhar-type equations 
as follow 

• 1 2 
Kii - t Yn 

f^i = ^ Yii Y 21 

" " O' ^21 

Y21 ^ K21 ) Yii- ^ Y21 



Fig. 5 shows the plot of the Kalman gain for this problem where we 
assumed that 

C - 0.01 Farad 

R - 10.0 Ohm 

L = 1.0 Henry 
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Case II ; (High Initial Uncertainty) 

If the initial state X(to) is not known, then we may assume that 
EEXjjXo] = Hq = « (2-40) 

In this case what we should do for both Riccati-type and 
Chandrasekhar-type equations is work with P’^t) instead of P(t). 
Usirfg equation (2-10) we can get the relation 

gf [P‘'(t)] = -P’'(t)Ht)P'^t) (2-41) 

Substitute P(t) into the above equation, then we have 


-1 [P*‘^(t)3 = -P"\t)F- - Fl^p“\t) + H'H - P*^t) 
dt c 

GQ^G'p"'(t) 

where we assumed that 

Fc = F - GCH Qc = Q - CC" 

And if we let 

S(t) = p"'(t) 

R(t) = P'\t)X(t) 

Then the Kalman-filter equation can be written as 
S(t) = P''(t)X(t) + g| [P'‘(t)]X(t) 


(2-42) 

(2-43) 

(2-44) 

(2-45) 

(2-46) 


R(t) = -IFj + GQ(.G'S(t)rR(t) + H'Y(t). R(t ) = 0 (2-47) 


0 


or 
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And from (2-42) we have 

§{t) = -S(t)Fc - FcS(t) - S(t)GQcG's(t) + h'h. 

s(t„) = p-*(to) (2-48) 

The equation (2-44) - (2-48) are called "information-filter" form of 

the Kalman filter. And whenever we have problems with large ini‘'ial 

uncertainty, it is convenient to use them because 
Hq = P(to) = 0 -i- S(to) = 0 

To obtain Chandrasekhar-type equation, let 

6Qj.G' = LjLj (2-49) 

Now let us define an axn matrix A(t) such that 

A(t) = l; S(t) (2-50) 

Using the above equation, then the Kalman-fil ter equation (2-47) will 
become 

R(t) = -IFj, + LiA(t)]"R(t) + H''Y(t), R(to) = 0 (2-51) 

and it can be shown that A(t) in (2-51) can be obtained via the fol- 

lowing equations 

A(t) = UB(t)B"(t), A(to) = 0 (2-5?) 

B(t) = -[Fc + LjA(t)r B(t), B(to) = h' (2-53) 
where S(t) can be calculated via the equation 

S(t) = B(t)B"(T)dT 
0 


(2-54) 
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Case III : (Stationary Processes) 

In this case we assume that 

Elxoxo] = n = n (2-55) 

Since F is a stability matrix and the eigenvalues of F have negative 
real parts, then as |t - to| -> «» Z(*) will reduce to a stationary 
process. 

It can be shown, as was first done by Doob (1944), that for Z(*) 
as in (2-3) - (2-6) with constant parameters 

ElZ(t) Z'(s)] = H EtX(t) X" (s)] h" (2-56) 

where we have the relation 


EtX(t) X"(s)] = “ ^^n(s). t > s 

^ n(t)e^ (s - t)^ t < s 


(2-57) 


Consider the relations 

^<l>(t,to) = F<j)(t,to) 

2^ 4* (i^ii-g) “ 4» ('t,tQ)F 

,g^ f(t.f)dT = f(t,t) + ^ f(t,x)dT 


(2-58) 


(2-59) 


(2-60) 


31 


where '}»(t.to) is the transition matrix, then it is not hard to show 
that the variance matrix n(t) - E[X{t)X'(t)] obeys a linear matrix 
differential equation 

i(t) = Fn(t) + n{t)p" + GQg", n(to) = Hq (2-61) 


and an explicit solution to the above differential equation is 


F(t-to) F'(t-to) . t F(t-T) F'(t-T) 

n(t) = e ® n e ° + / e ^ ' 


GQG e 


dr (2-62) 


to 


Now by previous assumption that F is a stability matrix, we see that 
as t goes to n(t) tends to a constant matrix if. In fact from (2-61) 

it is clear that if is the unique solution of the equation 


Fff + nr' + GQG" = 0 (2-63) 

This is called matrix Lyapunov equation. 

To obtain Chandrasekhar equation when Hq = F, we precede as 
fol 1 ow; 

D =^-(nH' + 6C)(nH" + 6C)" (2-64) 


and when a = P we can take 

Lj = Fh' + GC (2-65) 

so that the Chandrasekhar-type equations become 

k(t) = -Y2(t)Y2(t)H", K(to) ^ nH' + GC (2-66) 


* 

= IF - K(t)H]Yj(t), Yj(tj) = 


(24i7) 
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And also note that P(t) can be found via 

P(t) = -Y2(t)Y;(t). P(tp) = If (2-68) 

And therefore 

P(t) = if - /* YzCt) Y2 (t)(1t (2-69) 

0 

Case IV : (Stationary Processes with known Covariances) 

This case is actually very close to the usual assumptions of 
statistical communication theory. Let us assume that the stationary 
signal and noise processes are given not via a state-space model but 
by their covariance functions. So that we assume 

y(t) = Z(t) + v(t), t > to (2-70) 


where we have 

EIv(t)v'(s)] = I6{t - s) 
ElZ(t)v'‘(s)] ^0 s > t 


E[y(t)y'(s)] = I6(t - s) + K(t,s) 


where K(t,s) is given by 


K(t,s) - El 2(t)Z'(s) Z(t)v"{s) + v(t)Z-(s) 1 


Me 


F(t - s) 


N, t >■ s 


N 


"e^'ts “ t) M', t < s 


(2-71) 

(2-72) 

(2-73) 

(2-74) 
(2-7 5) 
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where: 

F is an (n x n) matrix 
M is an (p x n) matrix 
N is an (n x p) matrix 

To make this consistent with Case III we shall assume that 

N has rank a, a < m,n(n,p) (2-76) 

Note that if we assume that Z(*) and v(*) are uncorrelated, K{t,s) is 
a covariance function, but here we assume a (one-sided) dependence of 
signal on past signal and noise, therefore K(t,s) itself will not, in 
general, be a covariance we remark that such a one-sided dependence 
arises naturally in feedback control and communication problems. 

Generally it is not easy to find a solution for this problem with- 
out first having to determine models for Z(*) and y(*). Therefore let 


us assume the following equations: 


Z(t) = Ms{t) 

(2-77) 

e(t) = Fe(t) + K(t) [Y(t) - He(t)] 

(2-78) 

K(t) = N - nm' 

(2-79) 


where the nxn matrix i:(«) obeys the Riccati equation 


lit) = Fl(t) + I(t)F' + K(t)K'(t), ^(to) = 0 (2-80) 

Under the assumption that the matrices M, F, and N are constant, 
we may calculate K(-) via the following n(p + pj) Chandrasekhar-type 
equations: 
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K(t) = -Y(t)Y'(t)H' 

(2-81) 


Y(t) = (F - K(t)H) Y(t) 

(2-82) 


K(t„) = N Y(t„) = N, 

(2-83) 


where the initial conditions N and Nq satisfy 

NN" = n,n; 

(2-84) 


and 

Ng n X Pi matrix 

(2-85) 


Note that whenever a = p, where a is rank of 

NN', we have much 


simpler formula, since we can take 

Y(to) = K(tj,) - N 

(2-86) 


Also note that there Chandrasekhar equations have 

the same form as in 



Case I. 
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CHAPTER III 

CONNECTIOH BETWEEN THE CHANDP>ASEKHAR X, Y 
FUNCTIONS AND THE CHAT^DRASEKHAR-TYPE EQUATIONS 

The algorithm is now in a form where we can point out a close 
relationship to some famous equations obtained in astrophysics. It 
is shov/n that the Kalman gain function was intimately related to the 
basic X and Y functions of Chandrasekhar (1). 

Consider the stochastic process 

y(t) - Z(t) + v(t), 0 < t < T < “ (3-1) 

where y is a white Gaussian noise process with unit spectral intensity, 
•such that 

E[v(t)] = 0. E[v(t)v'(s)] = 6{t - s) (3-2) 

and Z(t) is such that 

ElZ(t)v'(s)] ^0. t < s (3-3) 

And V and Z are n-dimensional vectors. Also for simplicity we assume 
that 

ElZ{t)] ==^0, 0 < t < T 

Therefore the covariance K of the signal process Z is the nxn 
matrix 

K(t,s) = E[Z(t)Z'(s)I, 0 < t. s 5 T. (3-4) 

where we impose the condition that K(t,s) is continuous in t and s. 

A 

Let us define the quantity Z(t) to be the linear least squares estimate 
of Z(t) to be given the observation {y(s), 0 < s < t), 0 < t < T. 



36 


I I 


It Is now a well-known fact that 2{t) may be written as a linear func- 
tional of Y, 

2(t) = h{t,s)y(s)ds* (3-5) 

0 

where the weighting function h(t,s) satisfies the integral equation 

h(t,s) - K(t,s) - K(T,s)h(t,T)dT, 0 < s < t < T (3-6) 

0 

In view of (3-5), all the information necessary to characterize the 
optimal estimate is contained in the function h(t,s). Let us now 
turn our attention toward development of a procedure for obtaining h 
as a function of the observation time t. 

Consider the integral equation 

h(t,s) = K(t,s) - K(T,s)h(t,x)dT, 0 < s < t <T (3-7) 

0 

where the kernel K is subject to the previous assumptions. We now 
wish to impose the additional assumption that K may be represented in 
the form: 

K(t,s) = K(|t-sD « exp.(-jt-sja)w(a)da (3-8) 

0 

and, by suitable choice of w(sums of delta functions, for example), a 
broad class of important problems may be simply handled. To avoid 
trivial details, we shall assume to is a scalar function. 

We can show that an initial value representation for h is given 
in terms of the functions X, Y and J as 
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-Y(t,a) /* Y(t.a'Ma') d»' (3-9) 

0 

alllisi = -aYtt.c) -X(t,a) Y(t,a')^a')da' (3-10) 

oX 0 

aJ(s.t.g) , .cJ( 5 ,t,a)-X(t.ci) /’ J(s,t,o')u(o')da' (3-11) 

3t 0 

for 0<S<t<T, 0<a<1. 

x(0,a) = 1 (3-12) 

Y(0,a) - 1 (3-13) 

, 0(S,S,a) = X(S,a) (3-14) 

The function h is then given by 

h(t,s) == J(s,t,a )w(a )(Ja (3-15) 

n 


The two nonlinear differential equations (3-9) and (3-10 are 
generally known as Chandrasekhar's X and Y functions. 

Although knowledge of the function h(t,s) is sufficient to deter- 
mine the optimal estimate 2(t) via (3-5), it is computationally desir- 

A 

able to obtain an alternative form for Z involving a quadrature over 
a fixed interval, rather than the variable interval of (3-5). Such a 

A 

representation will also simplify the real-time calculation of Z. 

Z(t) = /^h(t,s)Y(s)ds. t > 0 
0 


Recall (3-5); 
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Making use of (3-15), we see that 

Z(t) = 0(s,t,a )w{a )do( y(s)ds (3-16) 

0 0 

Interchanging the order of integration in (3-16) gives 

Z(tj = i3{s,t,a )y(s)ds ai(a )da (3-17) 

0 0 

Introduce the new function L(t,a) as 

L(t,a) = A J(s,t,a)y(s)ds, t > 0, 0 < a £ 1 (3-18) 

0 

We can show that L actually satisfies the Cauchy problem: 

L(t,a) = -aL(t,a) +X(t,a){y(t) -I L(t,ot )o)(a ) 

9t 0 

da'} (3-19) 

L(t,0) = 0 (3-20) 

/V 

By using these results, the optimal estimate Z is then given by 

Z(t) - S L(t,a )cA)(a )da (3-21) 

0 

At this point it seems essential to point out that unless the 
function h(t,s) is desired for some reason other than calculating Z, 
we need not calculate it at all. In fact, (3-21) and (3-19) show 
that knowledge of the function X and L suffice to determine Z. This 
observation is of considerable computational importance, as well as 
of some analytical interest. 
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With some effort, we should be able to see that these equations 
are essentially the same as (2-77), (2-78), (2-81), {2-82), (2-86) if 
we make the assumption 

n 

W(o) = ^ a^fi(a - a^),ct^ > 0, -F = diag{a ^ . ,a^) 

1 

This is why the equation for K(*) and Y(*) are said to be of 
Chandrasekhar type. 


Further Examine of Case IV 

To make some connections with the radiative-transfer problem in 
which Chandrasekhar originally introduced the X(-) and Y(0 functions, 
we .shall examine Special Case Ilia further. Let us assume also that 
N has rank P, so that the relevant equations are (2-81), (2-82), and 
(2-86). Not let us define an impulse response function h(t,T) as 

Z(t) = /■ h(t.T)y(t)dT (3-22) 

0 

Then we can show that 


h(t,T) 


f MV>{t,T)K(T), 

0 . 


t > T 

t < T 


(3-23) 


where i/'(.,.) is the state-transition matrix of F-K(-)M, that is 


(F - K(t) M) ♦ (t.T). Mt.t) = I (3-24) 


Also note that h(.,,} satisfies the Wiener-Hopf type of equation 
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h(t.x) + h(t.s)K(s.TMs = K(t,x).t > T (3-25) 

0 

Now it is clear that from (3-23) and (2-7), (2-82) we can obtain the 
following relations 

h(t,t) = MK(t) (3-26) 

h{t,0) = Mi'(t,0)K(0) = MY(t) (3-27) 

Then we see that the equation (2-81) 

k(t) = -Y(t)Y'(t)M', is the same as 

- h(t,0) h"(t,0) (3-28) 


which is a result that follows from a resolvent identifying in the 
theory of integral equations. 

It is of interest to not that when the signal !(•) and noise 
v(*) are uncorrelated, i.e., when K(t,s) is a covariance function, we 
can identify 

MK(t) = E[Z (t) Z"(t)] (3-29) 

MY(t) = EIZ (t) 2'(0)1 (3-30) , 

where Z(t) = Z(t) - Z(t) is the instantaneous error in the estimate 
of Z(t). The relations (3-29) and (3-30) have been used to suggest 
filtering-theory analogs of some quantities that arise in radiative 
transfer, e.g., the reflection, transmission, and internal intensity 
functions; however, these analogies have to be treated with some care 
because the assumption that K(t.s) is a covariance function, which is 


necessary for (3-29) and (3-30), is not valid in the radiative- 
transfer problem. 
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CHAPTER IV 

PROOF OF THE CHANDRASEKHAR- TYPE FORMULAS 


The Chandrasekhar-type equations presented in previous chapters 
can be derived in several ways. We choose one that begins with 
knowledge of the Kalman filter solution. In this way, the manner in 
which the constancy of the model parameters can be exploited will per- 
haps standout most clearly. Also, all questions of existence and 
uniqueness can be resolved immediately by use of the known results for 
the Kalman Theory (5). 

For our convenience, let us restate the equations of the Kalman 
filter 

X(t) = FX(t) + K(t)[Y(t) - HX(t)] , X(to) = 0 (4-1) 


KU) * P(t) h' + GC 


(4-2) 


(4-3) 


P(t) « FP(t) + P(t)F" - K(t)K'(t) + GQG', 

P('to) “ ^0 

Let us define '^(t,to) to be the state-transition matrix of the 
matrix F - K(t)H. That is 4'(t*to) is the unique solution of the linear 


° IF - K(t)Hl Mt-td), = I (4-4) 

Before we proceed to derive our algorithm, let us consider the 


following lemma 


« 

Lemma: The derivative P(-) of the solution P{*} of the Riccati- 

type equation (2-12) can be written as follow 

^(t) = »I>(t,to} P (to) ^ (t.to) (4-5) 

where 

Htj) = D = FP(tj)+ P{t„)F"- K(t„)K'(tj) + GQG' (4-6) 

And if we let P(to) = Uo nnd K(to) = PoH* + GC we have: 

P(tj) = D = FHj + HjF' - (II(|K' + GC){n(,H' + GC)'+ GQG' (4-7) 
Proof: we attack (4-3) by differentiation to get 

P(t) *= FP(t; + P(t)F - K(t)i:'(t) - K(t)K"(t) (4-8) 

Also, differentiation of both sides of (4-2) yield: 

^t) = P(t) h" (4-9) 

Using equation (4-9), we can rearrange (4-8) as 

'p(t) « (F - K(t) H) P(t).+ P(t)(F - K(t) H)" (4-10) 

Now let us temporarily di regard the dependence of K(*) on P(*)» and 
just regard K(-) as some given function, then equation (4-10) is a 
linear homogeneous equation in P{*)» This equation can be solved to 
give 

P(t) - 'l^(t,to) P(t ) il^'(t.to) 

U can be easily checked by differentiation that the last expression 
for P(.) does indeed satisfy (4-10). 
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The basic idea underlying the rest of the arguments is that for 
certain choices of P(tg), P(tg), and therefore P{t) can have rank less 
than n. For example let us consider Case I. 

Proof of the Case I 

Since we assumed that Rg = 0, therefore (4- 7) becomes 
P(to) = 

where 

Qc « Q - CC' 

And the basic relation (4-5) becomes 

P(t) = if{t.t;,)G()(,GV.(t,t5) (4-n) 

Since GQ^G^ is nonegetive-definite, we can factor it as 

GQ^G^ LjLj , Lj = n X o matrix (4-12) 

where 

rank of GQ^.G^ = a 

Note that the dimension of Lj is unique. Now if we define 

Yl(t) = tKt.tijjLi (4-13) 

Then by differentiation we see that 

Yi(t) = tF - K(t)H] Yi(t), Yi(to) = Li (4-14) 

Also from (4-11) - (4-13) we see that 
P(t) = Yj(t) Vl(t) 
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Then (4-18) can be written as 


P(t) = Yj(t) y'(t) - YjCt) Y'(t) 

(4-20) 

And also from (4-19) 


Y^(t) - (F - Kit) H) Y^(t), Y,(t ) = 

(4-21) 

Then' we see that 


K(t) = ^t)H' = IYi(t)Y,'(t) - Y2(t)Y2(t)l h' 

(4-22) 

K(to) “ = iI|,h' + GC 

(4-23) 


Equations (4-21) - (4-23) are just the Chandrasekhar-type equations* 
and the proof of the general formula is completed. 

Proof of the Case IV 

In this case we have covariance Information rather than a state 
model. Let us work with the Riccati equation (2-80). 

i(t) = Fl(t) + I(t)F' + K(t)K'(t), I(tj) = 0 
K(t) = N - 

and thus for I (t) we can have 

i (t) = (F - K(t)M) lit) + I it) (F - K(t)H)' (4-Z4) 

toWj if we define i|)(t,to)to be the state- transition matrix of 
F - K(t) M, then from (4-24) we have 


I(t) = il»(t,to) I (to) 'J^'(t»to) 


(4-25) 
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and from (4-9} we have 

R(t) = P(t)H' = Yi(t) Yj{t)H' (4-15) 

K(to) = P(to)H" + GC = GC (4-16) 

Equations (4-14) - (4-16) are simultaneous Chandrasekhar-type equa- 
tions for Case I. 

Proof of the General Formulas 

The proof in the general case is almost as simple. Since for 
any value of no, D = P(to) as given in equation (4-7) is symmetric, 
therefore, it can be factored as 

D = P(to) = hl'i - (4-17) 

Let us assume that 

a = rank of D 

and rank of LiLi is the number, say &, where & is the number of 
positive eigenvalues of D and the rank of L 2 L 2 is a - 3 . 

Here we can always choose the dimensions of Lj and L 2 so that 

« n X e, " n x(a - e). 

By these assumptions, now we can write (4-5) as 

P(t) = ii»(t,to)LiLi4» (t.to) - i}»(t.to)L 2 L 2 ^ (t.to) (4-18) 
Now if we define 

Yi(t) = if»(t,tQ)Li , 


i = 1.2 


(4-19) 
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and from (2-80) we get 

lit ) « K(to) K“(to) = ui (4-26) 

Now if N has rank P* < min(n.p)* we can write 

NN^ = NqNq, Njj an n X p, matrix 
Therefore from (4-25) 

i(t) = f'(t.to)NoNS''''(t.to) = Y(t) Y'(t) (4-Z7) 

where we can take 

Y(t) = »(t,to)N|, 

and therefore 

Y(t) = (F - K(t) M) Y(t), Y(to) = No (4-28) 

Note that when Pj = P we can identify N and Kq . Equations (4-26) - 
(4-28) are Chandrasekhar-type equations. 

Also note that we can verify the relation 

lit) = n(t) - P(t) 

so that 

lit) = 5 - P(t) = -P(t) 

And if K(t,s) as in (2-75) arises from ? model of the form (2-1) - 
(2-6), we can identify 

M * H, N(t) = n(t) H" + GC 




CONCLUSION 


A new class of algorithm for linear least squares estimation 
problem is presented. Ths new algorithm does have some computational 
saving when the number of observation vectors are much less than the 
dimension of the state space. 
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APPENDIX I 

Fortran IV Program for Solution of the Chandrasekhar-Type Equations 

Here we present the computer program, using the fourth-order 
Runge-Kutta method, to solve the Chandrasekhar-type equations of the 
example 3 in Cuapter II, Case I. 

The purpose of the Runge-Kutta method is to obtain an approxi- 
mate solution of a system of first order ordinary differential 
equations. 
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call PLrcC(-uOt-i.oj 
call exit 

LM 
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